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Abstract 

Motivated by recent interest for multi-agent systems and smart power grid architectures, we discuss 
the synchronization problem for the network-reduced model of a power system with non-trivial transfer 
conductances. Our key insight is to exploit the relationship between the power network model and 
a first-order model of coupled oscillators. Assuming overdamped generators (possibly due to local 
excitation controllers), a singular perturbation analysis shows the equivalence between the classic swing 
equations and a non-uniform Kuramoto model. Here, non-uniform Kuramoto oscillators are characterized 
by multiple time constants, non-homogeneous coupling, and non-uniform phase shifts. Extending methods 
from transient stability, synchronization theory, and consensus protocols, we establish sufficient conditions 
for synchronization of non-uniform Kuramoto oscillators. These conditions reduce to and improve upon 
previously-available tests for the standard Kuramoto model. Combining our singular perturbation and 
Kuramoto analyses, we derive concise and purely algebraic conditions that relate synchronization and 
transient stability of a power network to the underlying system parameters and initial conditions. 

I. Introduction 

The vast North American interconnected power grid is often referred to as the largest and most complex 
machine engineered by humankind. The various instabilities arising in such a large-scale power grid can 
be classified by their physical nature, the size of the uncertainty or disturbance causing the instability, 
or depending on the devices, processes, and the time necessary to determine the instability. All of these 
instabilities can lead and have led to blackouts of power grids [2], and their detection and rejection will 
be one of the major challenges faced by the future "smart power grid." 
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The envisioned future power generation will rely increasingly on renewable energy sources such as wind 
and solar power. Since these renewable power sources are highly stochastic, there will be an increasing 
number of transient disturbances acting on a power grid that is expected to be even more complex and 
decentralized than the current one. Thus, an important form of power network stability is the so-called 
transient stability [3], which is the ability of a power system to remain in synchronism when subjected 
to large transient disturbances. These disturbances may include faults on transmission elements or loss 
of load, loss of generation, or loss of system components. For example, a recent major blackout in Italy 
in 2003 was caused by tripping of a tie-line and resulted in a cascade of events leading to the loss of 
synchronism of the Italian power grid with the rest of Europe [2] . The mechanism by which interconnected 
synchronous machines maintain synchronism is a balance of their mechanical power inputs and their 
electrical power outputs depending on the relative rotor angles among machines. In a classic setting the 
transient stability problem is posed as a special case of the more general synchronization problem, which 
is defined over a possibly longer time horizon, for rotor angles possibly drifting away from their nominal 
values, and for generators subject to local excitation controllers aiming to restore synchronism. In order 
to analyze the stability of a synchronous operating point of a power grid and to estimate its region of 
attraction, various sophisticated algorithms have been developed [4], [5], [6], [7], [8], [9]. Reviews and 
survey articles on transient stability analysis can be found in [10], [11], [12], [13]. Unfortunately, the 
existing methods can cope only with simplified models and do not provide simple formulas to check if a 
power system synchronizes for a given system state and parameters. In fact, an open problem, recognized 
by [14] and not resolved by classical analysis methods, is the quest for explicit and concise conditions for 
synchronization as a function of the topological, algebraic, and spectral graph properties of the network. 

The recent years have witnessed a burgeoning interest of the control community in cooperative control 
of autonomous agent systems. Recent surveys and monographs include [15], [16], [17]. One of the basic 
tasks in a multi-agent system is a consensus of the agents' states to a common value. This consensus 
problem has been subject to fundamental research [18], [19] as well as to applications in robotic coordi- 
nation, distributed sensing and computation, and various other fields including synchronization. In most 
articles treating consensus problems the agents obey single integrator dynamics, but the synchronization 
of interconnected power systems has often been envisioned as possible future application [20] . However, 
we are aware of only one article [21] that indeed applies consensus methods to a power network model. 

Another set of literature relevant to our investigation is the synchronization of coupled oscillators [22], 
in particular in the classic model introduced by Kuramoto [23]. The synchronization of coupled Kuramoto 
oscillators has been widely studied by the physics [24], [25], [26] and the dynamical systems communities 
[27], [28], [29]. This vast literature with numerous theoretical results and rich applications to various 
scientific areas is elegantly reviewed in [30], [31]. Recent works in the control community [18], [19], 
[32], [33] investigate the close relationship between Kuramoto oscillators and consensus networks. 

The three areas of power network synchronization, Kuramoto oscillators, and consensus protocols are 
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apparently closely related. Indeed, the similarity between the Kuramoto model and the power network 
models used in transient stability analysis is striking. Even though power networks have often been 
referred to as systems of coupled oscillators, the similarity to a second-order Kuramoto-type model has 
been mentioned only very recently in the power networks community in [34], [35], [36], where only 
qualitative simulation studies for simplified models are carried out. In the coupled-oscillators literature, 
second-order Kuramoto models similar to power network models have been analyzed in simulations and 
in the continuum limit; see [31] and references therein. However, we are aware of only two articles 
referring to power networks as possible application [22], [37]. In short, the evident relationship between 
power network synchronization, Kuramoto oscillators, and consensus protocols has been recognized, but 
the gap between the first and the second two topics has not been bridged yet in a thorough analysis. 

There are three main contributions in the present paper. As a first contribution, we present a coupled- 
oscillator approach to the problem of synchronization and transient stability in power networks. Via a 
singular perturbation analysis, we show that the transient stabiUty analysis for the classic swing equations 
with overdamped generators reduces, on a long time-scale, to the problem of synchronizing non-uniform 
Kuramoto oscillators with multiple time constants, non-homogeneous coupling, and non-uniform phase- 
shifts. This reduction to a non-uniform Kuramoto model is arguably the missing link connecting transient 
stability analysis to networked control, a link that was hinted at in [14], [20], [34], [35], [36], [22], [21]. 

Second, we give novel, simple, and purely algebraic conditions that are sufficient for synchronization 
and transient stability of a power network. To the best of our knowledge these conditions are the first 
ones to relate synchronization and performance of a power network directly to the underlying network 
parameters and initial state. Our conditions are based on different and possibly less restrictive assumptions 
than those obtained by classic analysis methods [4], [5], [6], [7], [8], [9], [10]. We consider a network- 
reduced model of a power network, and do not make any of the following common or classic assumptions: 
we do not require the swing equations to be formulated in relative coordinates accompanied by a uniform 
damping assumption, we do not require the existence of an infinite bus, and we do not require the transfer 
conductances to be "sufficiently small" or even negligible. On the other hand, our results are based on the 
assumption that each generator is strongly overdamped, possibly due to internal excitation control. This 
assumption allows us to perform a singular perturbation analysis and study a dimension-reduced system. 
Due to topological equivalence, our synchronization conditions hold locally even if generators are not 
overdamped, and in the application to real power networks the approximation via the dimension-reduced 
system is theoretically well-studied and also applied in the power industry [38]. Our synchronization 
conditions are based on an analytic approach whereas classic analysis methods [4], [5], [7], [8], [9], [10] 
rely on numerical procedures to approximate the region of attraction of an equilibrium by level sets of 
energy functions and stable manifolds. Compared to classic analysis methods, our analysis does not aim 
at providing best estimates of the region of attraction or the critical clearing time. Rather, we approach 
the open problem [14] of relating synchronization to the underlying network structure. For this problem. 
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we derive sufficient and purely algebraic conditions that can be interpreted as "the network connectivity 
has to dominate the network's non-uniformity, the network's losses, and the lack of phase cohesiveness." 

Third and final, we perform a synchronization analysis of non-uniform Kuramoto oscillators, as an 
interesting mathematical problem in its own right. Our analysis combines and extends methods from 
consensus protocols and synchronization theory. As an outcome, purely algebraic conditions on the 
network parameters and the system state establish the phase cohesiveness, frequency synchronization, and 
phase synchronization of the non-uniform Kuramoto oscillators. We emphasize that our results do not 
hold only for non-uniform network parameters but also in the case when the underlying coupling topology 
is not a complete graph. When our results are specialized to classic (uniform) Kuramoto oscillators, they 
reduce to and even improve upon various well-known conditions in the literature on the Kuramoto model 
[25], [24], [19], [32], [33], [39]. In the end, these conditions guaranteeing synchronization of non-uniform 
Kuramoto oscillators also suffice for the transient stability of the power network. 

Paper Organization: This article is organized as follows. The remainder of this section introduces 
some notation, recalls preliminaries on algebraic graph theory and differential geometry, and reviews 
the consensus protocol and the Kuramoto model of coupled oscillators. Section II reviews the problem 
of transient stability analysis. Section III introduces the non-uniform Kuramoto model and presents the 
main result of this article. Section IV translates the power network model to the non-uniform Kuramoto 
model whose synchronization analysis is presented in Section V. Section VI provides simulation studies 
to illustrate the analytical results. Finally, some conclusions are drawn in Section VII. The appendix in 
Section VIII contains different synchronization conditions and estimates for the phase cohesiveness that 
can be derived alternatively to the ones presented in Section V. 

Vector and matrix notation: Given an n-tuple (xi, . . . , rc„), diag(xj) G M"^" is the associated diagonal 
matrix, x G M" is the associated column vector, Xmax and x^m are the maximum and minimum elements, 
and ||x||2 and ||x||^ are the 2- and oo-norm. Let 1„ and 0^ be the vectors of I's and O's of dimension 
n. Given two non-zero vectors x G and y G M", the angle Z{x,y) G [0,7r/2] between them satisfies 
cos(Z(x,y)) = x-^y/dl^ll ||y||). Given an array {Aij} with i,j G {1, . . . ,n}, A G M"^" is the associated 
matrix with ylmax = niaxj j{j4jj} and Amin = m.mij{Aij}. Given a total order relation among the indices 
let diag{Aij) denote the corresponding diagonal matrix. 

Graph theory: A weighted directed graph is a triple Q = (V, £, A), where V = {1, . . . , n} is the set of 
nodes, C V x V is the set of directed edges, and A G R"^" is the adjacency matrix. The entries of A 
satisfy aij > for each directed edge G £ and are zero otherwise. Any nonnegative matrix A induces 
a weighted directed graph Q. The Laplacian of Q is the n x n matrix L{aij) := diag(^"^^ Uij) — A. In 
the following, we assume that A = A^, that is, Q is undirected. In this case, the graph Q is fully described 
by the elements aij with i > j. If a number k G {1,...,|<S|} and a weight Wk = aij is assigned to any 
of these edges with i > j, then the incidence matrix H G rI^I^" is defined component-wise as 

Hki = 1 if node / is the sink node of edge k and as H^i = — 1 if node I is the source node of edge k; all 
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Other elements are zero. The Laplacian equals then the symmetric matrix L{aij) = H'^ diag{wk)H. If Q 
is connected, then ker(//) = ker(L(ajj)) = span(l„), all n — 1 remaining non-zero eigenvalues of L{Q) 
are strictly positive, and the second-smallest eigenvalue X2{L{aij)) is called the algebraic connectivity of 
Q and, for a complete and uniformly weighted graph (ajj = 1 for all i ^ j), it satisfies X2{L{aij)) = n. 

Geometry on the n-torus: The torus is the set = ] — vr, +7r], where — vr and +7r are associated with 
each other, an angle is an element ^ G T^, and an arc is a connected subset of T^. The product set 
T" is the n-dimensional torus. With slight abuse of notation, let \9i — ^2! denote the geodesic distance 
between two angles 9i € and 6*2 G T^. For 7 G [0,7r], let A(7) C T" be the set of angle arrays 
(6*1, . . . , 6'„) G T" such that there exists an arc of length 7 containing all 6*1, . . . , 6*^ in its interior. Thus, 
an array of angles 9 G A(7) satisfies maxj jgji^ „| \9i — 9j\ < 7. For 7 G [0,7r], we also define A(7) 
to be the union of the set {6 G T" | 9i = 9j, i,j G {1, • • • ,n}} and the closure of the open set A (7). 

For a rigorous definition of the difference between angles (i.e., points on the torus), we restrict our 
attention to angles contained in an open half-circle: for angles 9i, 92 with |0i — 6*2] < vr, the difference 
9i — 92 is the number in ] — vr, 7r[ with magnitude equal to the geodesic distance |0i — ^2! and with positive 
sign iff the counter-clockwise path length connecting 61 and 92 on is smaller than the clockwise path 
length. Finally, we define the multivariable sine sin : T" — )• [0, 1]" by sin(j;) = (sin(xi), . . . ,sin(xn)) 
and the sine function sine : M — )• M by sinc(x) = sin(x)/x. 

Review of the Consensus Protocol and the Kuramoto Model: In a system of n autonomous agents, 
each characterized by a state variable xi G M, one of the most basic tasks is to achieve a consensus on 
a common state value, that is, all agent states Xi{t) converge to a common value Xqo G M as t — )• 00. 
Given a graph Q with adjacency matrix A describing the interaction between agents, a simple, linear, 
and continuous time algorithm to achieve consensus on the agents' state is the consensus protocol 



In vector notation the consensus protocol (1) takes the form x = —L{aij)x, which directly reveals the 
dependence of the consensus protocol to the underlying graph Q. 

A well-known and widely used model for the synchronization among coupled oscillators is the Ku- 
ramoto model, which considers n coupled oscillators with state 6i G with the dynamics 



where K is the coupling strength and is the natural frequency of oscillator i. Unlike for the consensus 
protocol (1), different levels of consensus or synchronization can be distinguished for the Kuramoto 
model (2): The case when all angles 9i{t) converge to a common angle 6*00 G as t — )• 00 is referred 
to as phase synchronization and can only occur if all natural frequencies are identical. If the natural 
frequencies are non-identical, then each phase difference 9i{t) — 9j{t) can converge to a constant value, 
but this value is not necessarily zero. A solution 9 : M>o — )■ T" to the Kuramoto model (2) is phase 
cohesive if there exists a length 7 G [0, it[ such that 9{t) G A(7) for all t > 0, i.e., at each time t there 
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exists an arc of length 7 containing all angles 6i{t). A solution 6 : M>o — )• T" achieves exponential 
frequency synchronization if all frequencies 6i{t) converge exponentially fast to a common frequency 
^00 £ as t — )• 00. Finally, a solution 9 : ]R>o — )■ T" achieves exponential synchronization if it is phase 
cohesive and it achieves exponential frequency synchronization. In this case, all phases become constant 
in a rotating coordinate frame with frequency ^00 > and hence the terminology phase locking is sometimes 
also used in the literature. 

II. Models and Problem Setup in Synchronization and Transient Stability Analysis 

A. The Mathematical Model of a Power Network 

In a power network with n generators we associate to each generator its internal voltage Ei > 0, its 
active power output Pe,i> its mechanical power input P^^i > 0, its inertia Mj > 0, its damping constant 
Di > 0, and its rotor angle 9i measured with respect to a rotating frame with frequency /q. All parameters 
are given in per unit system, except for Mj and Di which are given in seconds, and /o is typically given 
as 50 Hz or 60 Hz. The rotor dynamics of generator i are then given by the classic constant-voltage 
behind reactance model of interconnected swing equations [11], [40], [41] 

MiOi = Pm,i - EfGu - DA - Pe,i, i G {1, . . . , n} . 

Under the common assumption that the loads are modeled as passive admittances, all passive nodes of a 
power network can be eliminated (c.f. Kron reduction [42]) resulting in the reduced (transfer) admittance 
matrix Y = 1"-^ G C"^", where Ya is the self- admittance of generator i and '^{Yij) > and Q{Yij) > 0, 
i j, are the transfer conductance and (inductive) transfer susceptance between generator i and j in per 
unit values. With the power-angle relationship, the active output power Pe,i is then 

Given the transfer admittance Yij between generator i and j, define the magnitude |yij| > and the phase 
shift ifij = arctan(K(yij)/9(lij)) G [0,7r/2[ depicting the energy loss due to the transfer conductance 
di{Yij). Recall that a lossless network is characterized by zero phase shifts. Furthermore, we define 
the natural frequency uji := Pm,i — Ef^{Yii) (effective power input to generator i) and the coupling 
weights Pij := EiEj\Yij\ (maximum power transferred between generators i and j) with Pa := for 
i G {1, . . . , n}. The network-reduced power system model can then be formulated compactly as 

En 
,^^Pijsm{ei-ej + ipij). (3) 

Typically, a dynamical model for the internal voltage of generator i is given as Ei = Ei{Ei, Ui, 6i — Oj), 
where n, is the field excitation and can be used as a control input [43]. Higher order electrical and flux 
dynamics can be reduced [44] into an augmented damping constant Di in equation (3). The generator's 
internal excitation control essentially increases the damping torque towards the net frequency and can also 
be reduced into the damping constant Di [40], [44]. It is commonly agreed that the classical model (3) 
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captures the power system dynamics sufficiently well during the first swing. Thus, we omit higher order 
dynamics and control effects and assume they are incorporated into the model (3). We remark that all our 
results are also valid if Ei = Ei{t) is a smooth, bounded, and strictly positive time-varying parameter. 
A frequency equilibrium of (3) is characterized by ^ = and by the (reduced) real power flow equations 

En 
Pij sin(0i - Bj + ipij) = 0, iG{l,...,n}. (4) 

depicting the power balance. More general, the generators are said to be in a synchronous equilibrium if 
all angular distances 1 9i — 9j \ are constant and bounded (phase cohesive) and all frequencies are identical 
9i = 9j. Exponential synchronization is then understood as defined before for the Kuramoto model (2). 

In order to analyze the synchronization problem, system (3) is usually formulated in relative coordinates 
[45]. To render the resulting dynamics self-contained, uniform damping is sometimes assumed, i.e., Di/Mi 
is constant. Some other times, the existence of an infinite bus (a stationary generator without dynamics) 
as reference is postulated [4], [10]. We remark that both of these assumptions are not physically justified 
but are mathematical simplifications to reduce the synchronization problem to a stability analysis. 

B. Review of Classic Transient Stability Analysis 

Classically, transient stability analysis deals with a special case of the synchronization problem, namely 
the stability of a post-fault frequency equihbrium, that is, a new equilibrium of (3) arising after a change 
in the network parameters or topology. To answer this question various sophisticated analytic and numeric 
methods have been developed [10], [11], [12], [13], which typically employ the Hamiltonian structure of 
system (3). Since in general a Hamiltonian function for model (3) with non-trivial network conductance 
'^{Yij) > (or equivalently (pij > 0) does not exist [46], early transient stability approaches neglect the 
phase shifts ipij [4], [6], [10]. In this case, the power network model (3) takes form 

M9 = -D9 -VU{9f , (5) 

where V is the gradient and U :] — vr, vr]" — R is the potential energy given up to an additive constant by 

f^(^) = - EL ('^^^^ + E"=i (1 - - ^i))) • (6) 

When system (5) is formulated in relative or reference coordinates (that feature equilibria), the energy 
function {9, 9) ^ (1/2) 9'^M9 + U{9) serves (locally) as a Lyapunov function. In combination with the 
invariance principle, we clearly have that the dynamics (5) converge to ^ = and the largest invariant 
zero level set of VC/(0). In order to estimate the region of attraction of a stable equilibrium, algorithms 
such as PEBS [4] or BCU [7] consider the associated dimension-reduced gradient flow 

9 = -VU{9)'^. (7) 

Then {9*, 0) is a hyperbolic type-A; equilibrium of (5), i.e., the Jacobian has k stable eigenvalues, if and 
only if 9* is a hyperbolic type— A; equilibrium of (7), and if a generic transversality condition holds, then 
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the regions of attractions of both equilibria are bounded by the stable manifolds of the same unstable 
equilibria [4, Theorems 6.2-6.3]. This topological equivalence between (5) and (7) can also be extended 
to "sufficiently small" transfer conductances [7, Theorem 5.7]. For further interesting relationships among 
the systems (5) and (7), we refer to [4], [7], [12], [10]. Other approaches to lossy power networks with 
non-zero transfer conductances compute numerical energy functions [5] or employ an extended invariance 
principle [9]. Based on these results computational methods were developed to approximate the stability 
boundaries of (5) and (7) by level sets of energy functions or stable manifolds of unstable equilibria. 

To summarize the shortcomings of the classical transient stability analysis methods, they consider 
simplified models formulated in reference or relative coordinates (with uniform damping assumption) 
and result mostly in numerical procedures rather than in concise and simple conditions. For lossy power 
networks the cited articles consider either special benchmark problems or networks with "sufficiently 
small" transfer conductances. To the best of our knowledge there are no results quantifying this smallness 
of K(lij) or ipij for arbitrary networks. Moreover, from a network perspective the existing methods do not 
result in explicit and concise conditions relating synchronization to the network's state, parameters, and 
topology. The following sections will address these questions quantitatively via purely algebraic tests. 

III. The Non-Uniform Kuramoto Model and Main Synchronization Result 

A. The Non-Uniform Kuramoto Model 

As we have already mentioned, there is a striking similarity between the power network model (3) 
and the Kuramoto model (2). To study this similarity, we define the non-uniform Kuramoto model by 



where we assume that the parameters satisfy Di > 0, G M, Pij > 0, and ipij £ [0,7r/2[, for all 
hj £ {1; • • • ) '^l' ^ 7^ j; by convention. Pa and ifu are set to zero. System (8) may be regarded as a 
generalization of the classic Kuramoto model (2) with multiple time-constants Di and non-homogeneous 
but symmetric coupling terms Pij and phase shifts ipij. The non-uniform Kuramoto model (8) will serve 
as a link between the power network model (3), the Kuramoto model (2), and the consensus protocol (1). 

Remark III.l (Second-order systems and their first-order approximations:) The non-uniform Ku- 
ramoto model (8) can be seen as a long-time approximation of the second order system (3) for a small 
"inertia over damping ratio" Mi/ Di. Note the analogy between the non-uniform Kuramoto model (8) and 
the dimension-reduced gradient system (7) studied in classic transient stability analysis to approximate the 
stability properties of the second-order system (5) [4], [7], [10]. Both models are of first order, have the 
same right-hand side, and differ only in the time constants Di. Thus, both models have the same equilibria 
with the same stability properties and with regions of attractions bounded by the same separatrices [4, 
Theorems 3.1-3.4]. The reduced system (7) is formulated as a gradient-system to study the stability of 
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the equilibria of (7) (possibly in relative coordinates). The non-uniform Kuramoto model (8), on the other 
hand, can be directly used to study synchronization and reveals the underlying network structure. □ 

B. Main Synchronization Result 

We can now state our main result on the power network model (3) and the non-uniform Kuramoto model (8). 

Theorem III.2 (Main synchronization result) Consider the power network model (3) and the non- 
uniform Kuramoto model (8). Assume that the minimal lossless coupling of any oscillator to the network 
is larger than a critical value, i.e., 

Tmin := n mill | ^ cos{Lpij) \ > Tcriticai := r ^ f™5^^ "FT ~ ^ + ^ max ^ sin(v3y) 

i^i I A J COs((/?max) ^ A Dj ie{l,...,n} ^ A 

(9) 

Accordingly, define 7min £ [0,7r/2 — ipmax[ cifid 7max G ]'^/2,7r] as unique solutions to the equations 

Sin(7min) = Sin(7max) = COs((/7max) rcritical/rmin- 

For the non-uniform Kuramoto model, 

1) phase cohesiveness.' the set A(7) is positively invariant for every 7 G [7min>7max]. and each 
trajectory starting in A(7max) reaches A(7,jiin); and 

2) frequency synchronization: /or every 0(0) G A(7niax)> the frequencies 9i{t) synchronize exponen- 
tially to some frequency 6^0 £ [6*111111(0), 0max(O)]. 

For the power network model, for all 0(0) G A(7max) cind initial frequencies Oi(0), 

3) approximation error: there exists a constant e* > such that, if e := Mmax/^min < e*. then the 
solution {9{t),9{t)) of (3) exists for all t > and it holds uniformly in t that 

{e^{t) - 9n{t)) = {9i{t) - 9n{t)) + 0(e), Vi > 0, i G {1, . . . , n - 1}, 

(10) 

9{t) = D-^Q{9{t)) + 0{e), Vt > , 

where 9{t) is the solution to the non-uniform Kuramoto model (8) with initial condition 9{0) = 9(0), 
and D^^Q{9) is the power fiow (4) scaled by the inverse damping D^^; and 

4) asymptotic approximation error: there exists e and (^max sufficiently small, such that the 0{e) 
approximation errors in equation (10) converge to zero as t ^ 00. 

The proof of Theorem III.2 is based on a singular perturbation analysis of the power network model (3) 
(see Section IV) and a synchronization analysis of the non-uniform Kuramoto model (8) (see Section V) 
and will be postponed to the end of Section V. We discuss the assumption that the perturbation parameter 
e needs to be small separately in the next subsection and state the following remarks to Theorem III.2: 

Remark III.3 (Physical interpretation of Theorem III.2:) The right-hand side of condition (9) states 
the worst-case non- uniformity in natural frequencies (the difference in effective power inputs at each 
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generator) and the worst-case lossy coupling of a generator to the network (Pij sin{ipij) = EiEj^(Yij) 
reflects the transfer conductance), both of which are scaled with the rates Di. The term cos((^max) = 
sin(7r/2 — ipmsix) corresponds to phase cohesiveness in A(7r/2 — (pmax), which is necessary for the latter 
consensus-type analysis. These negative effects have to be dominated by the left-hand side of (9), which 
is a lower bound for minjj^"^^ (Pjj cos(93jj)/L'j) }, the worst-case lossless coupling of a node to the 
network. The multiplicative gap rcriticai/rmin between the right- and the left-hand side in (9) can be 
understood as a robustness margin that additionally gives a practical stability result determining the 
admissible initial and the possible ultimate lack of phase cohesiveness in A(7min) and A(7niax)- 

In summary, the conditions of Theorem III.2 read as "the network connectivity has to dominate the 
network's non-uniformity, the network's losses, and the lack of phase cohesiveness." In Theorem III.2 
we present the scalar synchronization condition (9), the estimate for the region of attraction A(7max)> 
and the ultimate phase cohesive set A(7min). In the derivations leading to Theorem III.2 it is possible 
to trade off a tighter synchronization condition against a looser estimate of the region of attraction, or a 
single loose scalar condition against n(n — 1)/2 tight pairwise conditions. These tradeoffs are explored in 
the appendix of this document. Finally, we remark that the coupling weights Pij in condition (9) are not 
only the reduced power flows but reflect for uniform voltages Ei and phase shifts ipij also the effective 
resistance of the original (non-reduced) network topology [42]. Moreover, condition (9) indicates at which 
generator the damping torque has to be increased or decreased (via local power system stabilizers) in 
order to meet the sufficient synchronization conditions. 

The power network model (3) inherits the synchronization condition (9) in the relative coordinates 
9i — On and up to the approximation error (10) which is of order e and eventually vanishes for e and 
ifmax sufficiently small. The relative coordinates can be shown to be well-posed (see Section IV). The 
convergence of the power network model only from almost all initial conditions is a consequence of the 
existence of saddle points in the non-uniform Kuramoto model. □ 

Remark III.4 (Refinement of Theorem III.2:) Theorem III.2 can also be stated for two-norm bounds 
on the parameters involving the algebraic connectivity (see Theorem V.5). For a lossless network, explicit 
values for the synchronization frequency and the exponential synchronization rate as well as conditions 
for phase synchronization can be derived (see Theorems V. 1 and V.IO). When specialized to the classic 
Kuramoto model (2), the sufficient condition (9) is improves the results [24], [25], [32], [33], [39], and 
it can also shown to be a tight bound. We refer the reader to the detailed comments in Section V. □ 

C. Discussion of the Perturbation Assumption 

The assumption that each generator is strongly overdamped is captured by the smallness of the 
perturbation parameter e = Mmax/^min- This choice of the perturbation parameter and the subsequent 
singular perturbation analysis (in Section IV) is similar to the analysis of Josephson arrays [26], coupled 
overdamped mechanical pendula [47], flocking models [48], and also classic transient stability analysis 
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[4, Theorem 5.2], [36]. In the Unear case, this analysis resembles the well-known overdamped harmonic 
oscillator, which features one slow and one fast eigenvalue. The overdamped harmonic oscillator exhibits 
two time-scales and the fast eigenvalue corresponding to the frequency damping can be neglected in the 
long-term phase dynamics. In the non-linear case these two distinct time-scales are captured by a singular 
perturbation analysis. In short, this reduction of a coupled-pendula system corresponds to the assumption 
that damping to a synchronization manifold and synchronization itself occur on separate time scales. 

In the application to realistic generator models one has to be careful under which operating conditions 
e is indeed a small physical quantity. Typically, Mi G [2s, 12s]/(27r/o) depending on the type of generator 
and the mechanical damping (including damper winding torques) is poor: Di G [1, 3]/(27r/o). However, 
for the synchronization problem also the generator's internal excitation control have to be considered, 
which increases the damping torque to Di G [10, 35]/(27r/o) depending on the system load [41], [40], 
[44]. In this case, e G 0(0.1) is indeed a small quantity and a singular perturbation approximation is 
accurate. In fact, the recent power systems literature discusses the need for sufficiently large damping to 
enhance transient stability, see [49], [50] and references therein. 

We note that simulation studies show an accurate approximation of the power network by the non- 
uniform Kuramoto model also for values of e G 0{\), i.e., they indicate that the threshold e* may be 
sizable. The theoretical reasoning is the topological equivalence discussed in Subsection II-B between 
the power network model (3) and the first-order model (7), which is again topologically equivalent to the 
non-uniform Kuramoto model (8), as discussed in Remark III. 1. The synchronization condition (9) on the 
non-uniform Kuramoto model (8) guarantees exponential stability of the non-uniform Kuramoto dynamics 
formulated in relative coordinates 9i — 9n, which again implies local exponential stability of the power 
network model (3) in relative coordinates. These arguments are elaborated in detail in the next section. 
Thus, from the viewpoint of topological equivalence. Theorem III.2 holds locally completely independent 
of e > 0, and the magnitude of e gives a bound on the approximation errors (10) during transients. 

The analogies between the power network model (3) and the reduced model (7), corresponding to the 
non-uniform Kuramoto model (8), are directly employed in the PEBS [4] and BCU algorithms [7]. These 
algorithms are not only scholastic but applied by the power industry [38], which additionally supports 
the validity of the approximation of the power network model by the non-uniform Kuramoto model. 
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IV. Singular Perturbation Analysis of Synchronization 
A. Time-Scale Separation of the Power Network Model 

In this section, we put the approximation of the power network model (3) by the non-uniform Kuramoto 
model (8) on solid mathematical ground via a singular perturbation analysis. The analysis by Tikhonov's 
method [51] requires a system evolving on Euclidean space and exponentially stable fixed points. In 
order to satisfy the assumptions of Tikhonov's theorem, we introduce two concepts. 

First, we introduce a smooth map from a suitable subset of T" to a compact subset of M"^^. For 
7 G [0,7r[, define the map grnd : A(7) — ;> Agi.nd(7) := G M""-*- | \5i\ < 7, maxjj \6i — Sj\ < 7, i,j e 
{1, . . . , n — 1}} that associates to the angles (^i, . . . , On) G ^(7) the array of angle differences 6 with 
components 6i = 6i — On, for i G {1, . . . ,n — 1}. This map is well defined, that is, 5 G Agrnd(7)> 
because = 1^1 — 6*^1 < 7 for all i G {1, ... ,n — 1} and |Ji — 5j| = \9i — 9j\ < 7 for all distinct 
i,j G {1, . . . ,n — 1}. Also, this map is smooth because 7 < vr implies that all angles take value in an 
open half-circle and their pairwise differences are smooth functions (see Section I). As a final remark, 
note that the angle differences ^1, . . . ,6n-i are well-known in the transient stability [7], [52] and in the 
Kuramoto literature [29], and we refer to them as grounded angles in the spirit of circuit theory. The 
sets A(7r) and Agmd(7r) as well as the map 0^5 = grnd(0) are illustrated in Figure 1. 




Fig. 1. Illustration of the map grnd : A (7) — > Agrnd(7)- Tlie map grnd can be tiiought of as as a symmetry -reducing projection 
from A (7) (illustrated as subset of to Ai,nid(7) (illustrated as subset of R^), where 6„ is projected to the origin 0. The set 
A(7) and the map grnd are invariant under translations on T" that is, under maps of the form (Si, ... , 9n) >->■ (0i+a, . . . S„ + a). 



Second, by formally computing the difference between the angles 9i and On, we define grounded 
Kuramoto model with state 5 G R"^^ by 



Di Dr. 



(11) 



June 28, 2011 



DRAFT 



13 



Lemma IV.l (Properties of grounded Kuramoto model) Let 7 G [0,7r[ and let 9 : M>o be a 

solution to the non-uniform Kuramoto model (8) satisfying 6(0) G ^(7). Let 6 : M>o — ?• M"^^ be the 
solution to the grounded Kuramoto model (11) with initial condition (5(0) = grnd(0(O)) G Agmd(7)- 
Then, 5{t) = grnd(9 (t)) for all t > 0, if any one of the two following equivalent conditions holds: 

1) phase cohesiveness: the angles 6{t) take value in A(7) for all time t > 0; and 

2) well-posedness: the grounded angles 5{t) take value in AgmdlT) for all time t > 0. 
Moreover, the following two statements are equivalent for any 7 G [0,7r[; 

3) exponential frequency synchronization: each trajectory of the non-uniform Kuramoto model 
satisfying the phase cohesiveness property 1) achieves exponential frequency synchronization; and 

4) exponential convergence to equilibria: each trajectory of the grounded Kuramoto model satisfying 
the well-posedness property 2) converges exponentially to an equilibrium point. 

Finally, every trajectory of the grounded Kuramoto model as in 4 ) satisfying the well-posedness property 
2) with 7 G [0,-7r/2 — fmax] converges to an exponentially stable equilibrium point. 

Proof: Note that, since both vector fields (8) and (11) are locally Lipschitz, existence and uniqueness 
of the corresponding solutions follows provided that the corresponding evolutions are compact. Now, 
assume that 1) holds, that is, 6{t) G A(7) (compact) for all t > 0. Therefore, 6{t) = grnd(0(f)) G 
Agrnd(7) for allt > 0. Also recall that the map grnd is smooth over A (7). These facts and the definition 
of the grounded Kuramoto model (11) imply that ^ grnd(0(t)) is well defined and identical to 6{t) for 
all t > 0. In turn, this implies that 6{t) = gvnd{9{t)) G Agrnd(7) holds for all positive times. 

Conversely, assume that 2) holds, that is, 6{t) G Agrnd(7) (compact) for all t > 0. Due to existence and 
uniqueness and since initially (5(0) = grnd(0(O)) with ^(0) G A(7), a set of angles 6{t) G A(7) can be 
associated to 6{t) G Agrnd(7) such that 6{t) = grnd(0(t)) for all t > 0. By construction of the grounded 
Kuramoto model (11), we have that 9{t) is identical to the solution to the non-uniform Kuramoto model 
(8) for all t > 0. Thus, statement 2) implies statement 1) and 6{t) = grnd(0(t)) for all t > 0. Having 
established the equivalence of 1) and 2), we do not further distinguish between 6{t) and grnd(0(t)). 

Assume that 3) holds, that is, all 6i{t) converge exponentially fast to some 600 G M. It follows 
that each 6i{t) = 9i{t) — 6n(t) converges exponentially fast to zero, and 6{t) = 6{0)+£6{T)dT converges 
exponentially fast to some 600 S Agrnd(7) due to property 2). Since the vector field (11) is continuous and 
lim (6{t),6{t)) = {600,0), the vector 600 is necessarily an equilibrium of (11), and property 4) follows. 

Assume that 4) holds, that is, all angular differences 6i{t) = 9i{t) — 6n{t) converge exponentially 
fast to constant values 5j 00 for i G {1, . . . , n — 1}. This fact and the continuity of the vector field in 
equation (11) imply that the array with entries (5j 00 is an equilibrium for (11) and that each frequency 
difference 6i{t) = 6i{t) — 9n{t) converges to zero. Moreover, because the vector field in equation (11) 
is analytic and the solution converges exponentially fast to an equiUbrium point, the right-hand side of 
equation (11) converges exponentially fast to zero and thus also the time-derivative of the solution, i.e., 
the array of all frequency differences, converges exponentially fast to zero. 
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To prove the final statement, assume that the non-uniform Kuramoto model (8) achieves frequency 
synchronization with synchronization frequency ^sync G and phase cohesiveness in A(7r/2 — tpmax)- 
Thus, when formulated in a rotating coordinate frame with zero synchronization frequency, all trajectories 
9i{t) — ^sync • t (mod 27r) necessarily converge to an equilibrium 6* G A(7r/2 — (pmax)- 

The Jacobian J {9*) of the non-uniform Kuramoto model is given by the Laplacian matrix with weights 
aij{9*) = {Pij / Di) cos{9* - 9* + (pij). For any 9* G A(7r/2 - (pm^y), the weights aij{9*) are strictly 
positive. In this case, it follows from the contraction property [53, Theorem 1] that the linearized dynamics 
9 = J {9*) ■ 9 converge from any initial condition in M" to a point in the diagonal vector space 1„. Hence, 
for any 9* G A(7r/2 — 9Jmax)> the matrix J{9*) has n — 1 stable eigenvalues and one zero eigenvalue with 
eigenspace 1„ corresponding to the translational invariance of the angular variable. 

Hence, any equilibrium manifold 9* G A(7r/2 — (^max) (of dimension one due to translational invariance) 
is exponentially stable w.r.t. to the non-uniform Kuramoto dynamics (8). The corresponding point 5* = 
grnd(6'*(t)) G Agrn(i(^/2 — '/'max) (the translational symmetry is removed) is an equilibrium of the 
grounded Kuramoto dynamics (11) (due to property 4)). Finally, since 9* is exponentially stable, it 
necessarily follows that 5* is an exponentially stable equilibrium point. ■ 

As mentioned in Remark III.l, system (8) may be seen a long-time approximation of (3), or spoken 
differently, it is the reduced system obtained by a singular perturbation analysis. A physically reasonable 
singular perturbation parameter is the worst-case choice of Mi/Di, that is, e := Mmax/^min- The 
dimension of e is in seconds, which makes sense since time still has to be normalized with respect 
to e. If we reformulate the power network model (3) in grounded angular coordinates with the state 
{5, 9) G R"^^ X M", then we obtain the following system in singular perturbation standard form 

^^6i=n{9):=9i-9n, i G {!,..., n-1}, (12) 

9i =gi{S, 9) := -Fi + ^ {^i - X^^-i ~ + '^^J^) ' ^ ^ {1, . . . , n} , (13) 

where Fi := {Di / Dmin) / {Mi / Mmax) and 5n '■= in equation (13). For e sufficiently small, the long-term 
dynamics of (12)-(13) can be approximated by the grounded Kuramoto model (11) and the power flow (4), 
where the approximation error is of order e and Fi determines its convergence rate in the fast time-scale. 

Theorem IV.2 (Singular Perturbation Approximation) Consider the power network model (3) written 
as the singular perturbation problem (12)-(13) with bounded initial conditions (5(0), ^(0)), and the 
grounded non-uniform Kuramoto model (11) with initial condition 5{^) and solution 5{t). Assume that 
there exists an exponentially stable fixed point 6oo of (11) and 6{0) is in a compact subset 0,s of its 
region of attraction. Then, for each Qs 

1) there exists £=„ > such that for all e < e*, the singular perturbation problem (12)-(13) has a 
unique solution {S{t, e),9{t, e)) for t > 0, and for all t > it holds uniformly in t that 

6{t,e) -6{t) = 0{e), and 9{t, e) - h{6{t)) - y{t/e) = 0{e) , (14) 
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where yi{t/e) := (^i(O) - /ii(5(0))) e"^'*/^ and hi{5) := Qi{6)/D,for i£{l,...,n}. 

2) For any ti, > 0, there exists e* < e^, such that for all t>tb and whenever e<e* it holds uniformly that 

e{t,e)-h{6{t)) = 0{e). (15) 

3) Additionally, there exist e and (/Jmax sufficiently small such that the approximation errors (14)-(15) 
converge exponentially to zero as t ^ oo. 

Proof: To prove statements 1) and 2) of Theorem IV.2 we will follow Tikhonov's theorem [51, 
Theorem 11.2] and show that the singularly perturbed system (12)-(13) satisfies all assumptions of [51, 
Theorem 11.2] when analyzing it on M"^^ x R" and after translating the arising fixed point to the origin. 

Exponential stability of the reduced system: The quasi-steady-state of (12)-(13) is obtained by solving 
gi{5, ^) = for 6, resulting in the unique (and thus isolated) root 9i = hi{6) = Qi{5)/ Di, i € {1, . . . , n}. 
The reduced system is obtained as 5i = fi{h{5)) = hi{5) — hn{6), i E {1, . . . , n — 1}, which is equivalent 
to the grounded non-uniform Kuramoto model (11). The reduced system is smooth, evolves on R"^^, 
and by assumption its solution 6{t) is bounded and converges exponentially to the stable fixed point 
6oo- Define the error coordinates Xi{t) = 6i{t) — 5i,oo, i G {l)---,'^ — 1} and the resulting system 
X = f{h{x + (5oo)) with state in R""^ and initial condition x{Q) = 5(0) — 5oo- By assumption, the 
solution x{t) is bounded and converges exponentially to the stable fixed point at x = 0. 

Exponential stability of the boundary layer system: Consider the error coordinate Hi = 9i — hi[5), 
which shifts the error made by the quasi-stationarity assumption 6i{t) ^ hi{6{t)) to the origin. After 
stretching time to the dimensionless variable t = t/e, the quasi-steady-state error is 

^ 2/. = 9i{5, y + h{5)) - f{y + h{5)) = -F, y, - My + h{6)) (16) 

with yi{0) = 9i{0) — hi{6{0)). By setting e = 0, (16) reduces to the boundary layer model 

^yi = -Fyi, y,{0) = 0,{O) - h,{6{0)) . (17) 

The boundary layer model (17) is globally exponentially stable with solution yi{t/e) = yi(0)e~^'*/'^, 
where yi{0) is in a compact subset of the region of attraction of (17) due to boundedness of (5(0), ^(0)). 

In summary, the singularly perturbed system (12)-(13) is smooth on R"^^ x R", and the origins of 
the reduced system (in error coordinates) x = f{h{x + 5oo)) and the boundary layer model (17) are 
exponentially stable (where Lyapunov functions are readily existent by converse arguments [5 1 , Theorem 
4.14]). Thus, all assumptions of [51, Theorem 11.2] are satisfied and statements 1) and 2) follow. 

To prove statement 3), note that 5{t) converges to an exponentially stable equilibrium point 6oo, and 
{6{t, e),9{t, e)) converges to an 0{e) neighborhood of ((5oo, ^(^oo))> where h{6oo) = 0. We now invoke 
classic topological equivalence arguments from the transient stability literature [4], [7]. Both the second 
order system (12)-(13) as well as the reduced system 5 = f{h{6)) correspond to the perturbed Hamiltonian 
system (8)-(9) in [7] and the perturbed gradient system (10) in [7], where the latter is considered in [7] 
with Di = I for all i. Consider for a moment the case when all ipij = 0. In this case, the reduced system 
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has a locally exponentially stable fixed point 5qo (for any Di > due to [4, Theorem 3.1]), and by [7, 
Theorem 5.1] we conclude that {6oo,0) is also a locally exponentially stable fixed point of the second 
order system (12)-(13). Furthermore, due to structural stability [7, Theorem 5.7, Rl], this conclusion 
holds also for sufficiently small phase shifts ipij. Thus, for sufficiently small e and (/7max> the solution of 
(12)-(13) converges exponentially to (5oo,0). Consequently, the approximation errors 6{t,e) — 6{t) and 
0{t,e) — h{S) as well as the boundary layer error y{t/e) converge exponentially to zero. ■ 
Theorem IV.2 can be interpreted geometrically as follows. The frequency dynamics of system (3) happen 
on a fast time-scale and converge exponentially to a slow manifold which can be approximated to first 
order by the scaled power flow D~^Q{6). On this slow manifold the long-term phase synchronization 
dynamics of system (3) are given by the non-uniform Kuramoto model (8). 

V. Synchronization of Non-uniform Kuramoto Oscillators 

This section combines and extends methods from the consensus and Kuramoto literature to analyze 
the non-uniform Kuramoto model (8). The role of the time constants Di and the phase shifts ^pij is 
immediately revealed when dividing by Di both hand sides of (8) and expanding the right-hand side as 

ei = ^- ^ i-^cos{ipij)sm{9i-Bj) + -^sm{ipij)cos{ei-9j)j. (18) 

The difficulties in the analysis of system (8) are the phase shift-induced lossy coupling {PijjDi) sin(99jj) 
X cos{6i — 9j) inhibiting synchronization and the non-symmetric coupling between an oscillator pair 
{i,j} via Pij/Di on the one hand and Pij/Dj on the other. Since the non-uniform Kuramoto model 
(8) is derived from the power network model (3), the underlying graph induced by P is complete and 
symmetric, i.e., except for the diagonal entries, the matrix P is fully populated and symmetric. For the 
sake of generality, this section considers the non-uniform Kuramoto model (8) under the assumption that 
the graph induced by P is neither complete nor symmetric, that is, some Pij are zero and P / P^ . 

A. Frequency Synchronization of Phase-Cohesive Oscillators 

Under the assumption of cohesive phases, the classic Kuramoto model (2) achieves frequency synchro- 
nization [32, Theorem 3.1], [39, Corollary 11]. An analogous result guarantees frequency synchronization 
of non-uniform Kuramoto oscillators (8) whenever the graph induced by P has a globally reachable node. 

Theorem V.l (Frequency synchronization) Consider the non-uniform Kuramoto model (8) where the 
graph induced by P has a globally reachable node. Assume that there exists 7 G [0, 7r/2 — (/?n,ax[ such that 
the (non-empty) set of bounded phase differences A(7) is positively invariant. Then for every 0(0) G ^(7), 

1) the frequencies 6i{t) synchronize exponentially to some frequency 600 G [0min(O), 0max(O)]; and 

2) if '.p max = and P = P^ , then 600 = ^ ■= ^i/ X]j '^^^ exponential synchronization rate 
is no worse than 

Afe =-\2{L{Pij)) cos(7) cos(Z(L>l, l)f/D^^x. (19) 
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In the definition of the convergence rate Afe in (19), the factor X2{L{Pij)) is the algebraic connectivity 
of the graph induced by P = , the factor 1/-Dmax is the slowest time constant of the non-uniform 
Kuramoto system (8), the proportionality Afe ~ 003(7) reflects the phase cohesiveness in A(7), and 
the proportionality Afe ~ cos(Z(Dl, 1))^ reflects the fact that the error coordinate 6 — 01 is for non- 
uniform damping terms Di not orthogonal to the agreement vector $71. For non-zero phase shifts a 
small signal analysis of the non-uniform Kuramoto model (18) reveals that the natural frequency of 
each oscillator diminishes as — Xljyi ^ij sin((/5jj). In this case, and for symmetric coupling P = P^ , 
the synchronization frequency 9 00 in statement 1) will be smaller than ^00 = in statement 2). When 
specialized to the classic Kuramoto model (2), statement 2) of Theorem V.l reduces to [32, Theorem 3.1]. 

Proof of Theorem V.l: By differentiating the non-uniform Kuramoto model (8) we obtain the 
following dynamical system describing the evolution of the frequencies 

DA = - V" , Pij cosiOi - Oj + ipij) {9, - Oj) . (20) 
at ^ — 'j=i 

Given the matrix P, consider a directed weighted graph Q induced by the matrix with elements Ojj = 
{Pij / Di) cos{6i — 6j + (pij). By assumption we have for every 9{0) G A(7) that 9{t) G ^(7) for all 
t > 0. Consequently the weights aij{t) = (Pij/Di) cos{9i{t) — 9j{t) + ipij) are non-degenerate, i.e., zero 
for Pij = and strictly positive otherwise for all t > 0. Note also that system (20) evolves on the tangent 
space of T", that is, the EucUdean space M". Therefore, the dynamics (20) can be analyzed as a linear 
time-varying consensus protocol for the velocities 9i with state-dependent Laplacian matrix L{aij): 

^J = -L{aij)9, (21) 

We analyze L{aij) as if it was a just time-varying Laplacian matrix L{aij{t)). At each time instance 
the matrix —L{aij{t)) is Metzler with zero row sums, and the weights aij{t) are bounded continuous 
functions of time that induce integrated over any non-zero time interval a graph with non-degenerate 
weights and a globally reachable node. It follows from the contraction property [53, Theorem 1] that 
9i{t) E [^min(O), ^max(O)] for alH > and 9i{t) converge exponentially to 9ao- This proves statement 1).^ 
In the case of zero shifts and symmetric coupling P = P^ the frequency dynamics (21) can be 
reformulated as a symmetric time- varying consensus protocol with multiple rates D as 

^^D9 = -L{w^j{t))9, (22) 

where L{wij{t)) is a symmetric time-varying Laplacian corresponding to a connected graph with strictly 
positive weights Wij{t) = Pij cos{9i — 9j). It follows from statement 1) that the oscillators synchronize ex- 
ponentially to some frequency ^oo- Since L{wij) is symmetric, it holds that 1^^ D9 = 0, or equivalently, 

'We remark that in the case of smoothly time-varying natural frequencies u)i{t) an additional term uj{t) appears on the right- 
hand side of the frequency consensus dynamics (21). If the natural frequencies are non-identical or not exponentially convergent 
to identical values, the oscillators clearly cannot achieve frequency synchronization and the proof of Theorem V.l fails. 
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DiOi{t) is a constant conserved quantity. If we apply this argument again at ^oo := linit-j>oo (^i{t), then 
we have Di6i{t) = DiOao, or equivalently, the frequencies synchronize exponentially to 9ao = ^■ 
In order to derive an explicit synchronization rate, consider the weighted disagreement vector 6 = 
d — Qln, as an error coordinate satisfying 1^-D(5 = lJ^D9 — iJ^D^ln = 0, that is, 6 lives in the weighted 
disagreement eigenspace of co-dimension 1„ and with normal vector D 1„. Since is constant and 
ker{L{wij)) = span(l„), the weighted disagreement dynamics are obtained from (44) in (^-coordinates as 

^^D6 = -L{w^j{t))6. (23) 

Consider the weighted disagreement function 5 ^ 6^D5 and its derivative along the dynamics (23) given by 

= -2 6'^L{wij{t))5. 

Since 6^ Din. = 0, it follows that 6 span(l„) and 6 can be uniquely decomposed into orthogonal 
components as (5 = {1^5 /n) 1„ + 5^, where 5±_ is the orthogonal projection of 5 on the subspace 
orthogonal to 1^. By the Courant-Fischer Theorem [54], the time derivative of the weighted disagreement 
function can be upper-bounded (point-wise in time) with the algebraic connectivity \2{L{Pij)) as follows: 

A s'^DS = -2 6lL{wij{t))5i_ = -{H6±f ■ diag(Pi, cos(^, - 0j)) ■ {Hd±) 

< -min|„}eHcos(0i-0,) : 9 G A{j)} ■ {HS^f -drngiPij) ■ {H6^) < -A2(L(Pi,)) 005(7) • . 

In the sequel, \\S±\\ will be bounded by \\6\\. In order to do so, let Ij^ = (1/ H^^x ||) 5_l be the unit 
vector that 6 is projected on (in the subspace orthogonal to 1„). The norm of 6± can be obtained as 
||<^±|| = ||<^"^1±|| = Pll cos(Z((5, 1_l)). The vectors 6 and 1± each live on (n — 1) -dimensional linear 
hyperplanes with normal vectors Din and !„, respectively, see Figure 2 for an illustration. The angle 
Z((5, 1±) is upper-bounded by max^ Z{6, 1±), which is said to be the dihedral angle and its sine is the 
gap between the two subspaces [54]. Since both hyperplanes are of co-dimension 1, we obtain the dihedral 
angle as the angle between the normal vectors Din and 1„, and it follows that Z{6, 1±) < Z{Dln, In) 
(with equality for n = 2). In summary, we have \\6\\ > \\6±\\ > \\6\\ cos(Z(L>l„, In))- 
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Finally, given Dmin ||^||^ < 6^D6 < Dmax and Afe as stated in equation (19), we obtain for the 
derivative of the disagreement function ^ 6^D6 < —2 \i^5^ D5. An application of the Bellman-Gronwall 
Lemma yields 6{tfD6{t) < 6{0)'^D6{0) e-2Aie(t) foj- t > 0. After reusing the bounds on S'^DS, we 
obtain that the disagreement vector 6{t) satisfies \\6{t)\\ < ^/D^^^^JD^\\5{0)\\ e~^'=(*) for all t>0. ■ 

B. Phase Cohesiveness 

The key assumption in Theorem V. 1 is that the angular distances are bounded in the set A(7r/2 — (^max)- 
This subsection provides two different approaches to deriving conditions for this phase cohesiveness 
assumption - the contraction property and ultimate boundedness arguments. The dynamical system 
describing the evolution of the phase differences for the non-uniform Kuramoto model (8) reads as 



Di Dj 



'^k=i -^k + fik) - sin(% -ek + fjk)^ , e {1, . . . ,n}. (24) 



Note that equation (24) cannot have a fixed point of the form 6i = 6j if the following condition is not met. 

Lemma V.2 (Necessary Condition on Synchronization) Consider the non-uniform Kuramoto model 
(8). For any two distinct i,j G {1, . . . ,n} there exists no solution of the form 9i{t) = 6j{t), t >0, if 



UJi UJj 



fe=i " ' ^ 

Condition (25) can be interpreted as "the coupling between oscillators i and j needs to dominate 
their non-uniformity" such that they can synchronize. For the classic Kuramoto model (2) condition (25) 
reduces to K < n/ (2(n — 1)) • [iOi — iOj), a necessary condition derived also in [32], [33], [25]. We remark 
that condition (25) is only a loose bound for synchronization since it does take into account the effect of 
lossy coupling induced by the phase shift ipij, which becomes obvious when expanding the sinusoidal 
coupling terms in (24) as in equation (18). Nevertheless, condition (25) indicates that the coupling needs 
to dominate the non-uniformity and possibly also disadvantageous effects of the lossy coupling. 

In order to show the phase cohesiveness 6{t) G A{tt/2 — (/?max)> the Kuramoto literature provides 
various methods such as quadratic Lyapunov functions [32], contraction mapping [33], geometric [24], 
or Hamiltonian arguments [27], [25]. Due to the non-symmetric coupling via the weights Pij/Di and 
the phase shifts ipij none of the mentioned methods appears to be easily applicable to the non-uniform 
Kuramoto model non-uniform Kuramoto model (8). A different approach from the literature on consensus 
protocols [18], [19], [39] is based on convexity and contraction and aims to show that the arc in which all 
phases are contained is of non-increasing length. A modification of this approach turns out to be applicable 
to non-uniform Kuramoto oscillators with a complete coupling graph and results in the following theorem. 

Theorem V.3 (Synchronization condition I) Consider the non-uniform Kuramoto-model (8), where the 
graph induced by P = P^ is complete. Assume that the minimal lossless coupling of any oscillator to 
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the network is larger than a critical value, i.e., 

( Pij , , 1 1 / iO'i iOj V ■\ Pij . , 

Tmin := n mm <^ — ^ cos((^jj) S > Tcriticai := t r max 7^ - 77- + 2 max > — ^ sm(v3jj) 

*7^i I A J COs{ipmax) ^ i^j Di Dj ie{l,...,n} ^ A 

(26) 

Accordingly, define 7min G [0,7r/2 — 93max[ 'Ji'^ 7max £ ]'^/2,vr] as unique solutions to the equations 

Sin(7min) = Sin(7max) = COs((/7max) rcritical/rmin- Then 

1) phase cohesiveness; the set A(7) is positively invariant for every 7 G [7min,7max]> and each 
trajectory starting in A(7max) reaches A(7niin); and 

2) frequency synchronization: /or every 6{0) £ the frequencies 9i{t) synchronize exponen- 
tially to some frequency 6^0 G [^min(O), 0max(O)]. 

Condition (26) is interpreted in Remark III.3. In essence, Theorem V.3 is based on the contraction 
property: the positive invariance of A(7) is equivalent to showing that all angles Oi{t) are contained 
in a rotating arc of non-increasing maximal length 7. This contraction analysis is similar to that of 
the consensus algorithms in [18], [19], [39], which derive their results on M". Throughout the proof of 
Theorem V.3 we comment on different possible branches leading to slightly different conditions than (26). 
These branches are explored in detail in the Appendix VIII. 

Remark V.4 (Reduction of Theorem V.3 to classic Kuramoto oscillators:) For the classic Kuramoto 
oscillators (2) the sufficient condition (26) of Theorem V.3 specializes to 

K > -ftTcritical •= '^max ~ '^min • (27) 

In other words, if > i^cnticah then for every 6'(0) G A(7max) the oscillators synchronize and are 
ultimately phase cohesive in A(7min), where 7max G ]7r/2, vr] and ^y^m G [0, 7r/2[ are the unique solutions 
to the equations sin(7niin) = sin(7max) = ^cnticai/^- To the best of our knowledge, condition (27) on the 
coupling gain K is the tightest explicit sufficient synchronization condition that has been presented in the 
Kuramoto literature so far. In fact, the bound (27) is close to the necessary condition for synchronization 
K > -ftTcriticai '^/(2(n — 1)) derived in Lemma V.2 and [32], [33], [25]. Obviously, for n = 2 condition (27) 
is necessary and sufficient for the onset of synchronization. Other sufficient bounds given in the literature 
scale asymptotically with n, e.g., [33, Theorem 2] or [32, proof of Theorem 4.1]. To compare our condition 
(27) with the bounds in [32], [24], [39], we note from the proof of Theorem V.3 that our condition can 
be equivalently stated as follows. The set A(7r/2 — 7), for 7 € ]0, it/2], is positively invariant if 

^ r^/ \ -^critical '^max '^min /ton 

K > K{'y) := = — . (28) 

cos(7j cos(7j 

Our bound (28) improves the bound K > K( j)n/2 derived in [32, proof of Theorem 4.1] via a quadratic 
Lyapunov function, the bound K > K(^)n/{n — 2) derived in [39] via contraction arguments similar to 
ours, and the bound derived geometrically in [24, proof of Proposition 1] that, after some manipulations, 
can be written in our notation as K > ^(7) cos((7r/2 — 7)/2)/ cos(7r/2 — 7). Our ongoing research also 
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reveals that the bound (27) is tight for a bimodal distribution of the natural frequencies € {wmax, '^min} 
and also satisfies the implicit consistency conditions in [29]. Thus, (27) is a necessary and sufficient con- 
dition for synchronization when the natural frequencies LOi are only known to be contained in [a;niin, Wmax]- 
We elaborate on this interesting circle of ideas in a separate publication [55]. 

In summary, condition (26) in Theorem V.3 improves the known [32], [33], [25], [24], [39] sufficient 
conditions for synchronization of classic Kuramoto oscillators, and it is a necessary and sufficient if the 
particular distribution of the natural frequencies cjj G [i^mim ^^max] is unknown. □ 

Proof of Theorem V3: We start by proving the positive invariance of A(7) for 7 G [0,7r]. Recall the 
geodesic distance between two angles on and define the non-smooth function [0, vr] by 

y(V') = maxllV'i - I i, j G {1, . . . ,n}}. 

By assumption, the angles 6i{t) belong to the set A(7) at time t = Q, that is, they are all contained 
in an arc of length 7 G [0,7r]. In this case, V{ip) can equivalently be written as maximum over a set 
of differentiable functions, that is, V{'4>) = maxjV'j — V'i I ^iJ S {!,..., n}}. The arc containing all 
angles has two boundary points: a counterclockwise maximum and a counterclockwise minimum. If we 
let Imax(^) (respectively /mm(V')) denote the set indices of the angles ipi, . . . ,iljn that are equal to the 
counterclockwise maximum (respectively the counterclockwise minimum), then we may write 

V{ll}) = iprn' - i^i', for all m' G /max(V') and / G lmin('0)- 

We aim to show that all angles remain in A(7) for all subsequent times t > 0. Note that 6{t) G A(7) 
if and only if V{9{t)) < 7. Therefore, A(7) is positively invariant if and only if V{9{t)) does not 
increase at any time t such that V{6{t)) = 7. The upper Dini derivative of V{6{t)) along the dynamical 
system (24) is given as in [19, Lemma 2.2] 

h\.0 a 

where m G Imax{G{t)) and £ G lrma{0{t)) are indices with the properties that 

Omit) = max{^„,(t) I m' G lmax(^(0)}, and ^^(t) = vcAn{e^.{t) \ 1' G /min(e(t))}. 
Written out in components (in the expanded form (18)) D^V{6{t)) takes the form 

D+V{e{t)) = ^ - ^ - V" {a^ksHOm{t) - Okit)) + a^ksmiOkit) - 9,{t))) 

Dm De ^k=l 

En 
{bmk cos(^^(t) - Okit)) - hk cos{9e{t) - Okit))) , (29) 
k=l 

where we used the abbreviations 0,^ := Pik cos{ipik) / Di and bik ■= Pik sm{ipik) / Di. The equality 
V{9{t)) = 7 implies that, measuring distances counterclockwise and modulo additional terms equal to 
multiples of 2-k, we have Omit) - 9i{t) = 7, < 6'm(i) - Okit) < 7, and < Okit) - Oeit) < 7. To 
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simplify the notation in the subsequent arguments, we do not aim at the tightest and least conservative 
bounding of the two sums on the right-hand side of (29) and continue as follows.^ 

Since both sinusoidal terms on the right-hand side of (29) are positive, they can be lower-bounded as 

amksin{9m(t)-0kit))+a(:ksm{9k{t)-0e{t)) > min {an,} (sm{em{t)-9k{t))+sm{0k{t)-0i{t))) 

ie{m,e}\{k} 

te{m,i}\{k} \ 2 J \ 2 

>2 min jc-t.) sinf — ) cos( — | = min ja,!.) sin(7) , 
- ie{m,e}\{k}^ ^ \2) \2) i6W}\{fc} 

where we applied the trigonometric identities sin(x)+sin(y) = 2 sin(^i^) cos(^^) and 2 sin(x) cos(y) = 
sin(x — y) + sin(x + y). The cosine terms in (29) can be lower bounded in A (7) as hmk cos{0m{t) — 
0k{t)) -hkcos{0t{t) - Okit)) > ~bmk ~ i>ik- In summary, D^V {0(t)) in (29) can be upper bounded by 

L>m iJl ^-^k=l i(z{m/}\{k} ^-^k=l ^-^k=l 



< max 



bJi UJj 



Pij 

n mm < — - 



1 V— \™ 

cos{ipij) > sin(7) + 2 max > bij , 
J ie{l,...,n}^ — 'i=i 



'J 

where we further maximized the coupling terms and the differences in natural frequencies over all possible 
pairs {m,£}. It follows that V{0{t)) is non-increasing for all 0{t) G A(7) if 

Tmin Sin(7) > COs(93max) Tcritical , (30) 

where Fmin and Fcnticai are defined in (26). The left-hand side of (30) is a strictly concave function of 7 E 
[0, vr]. Thus, there exists an open set of arc lengths 7 including 7* = vr/2 — (^^ax satisfying inequality (30) 
if and only if inequality (30) is true at 7* = vr/2 — (/Jmax with the strict inequality sign, which corresponds 
to condition (26) in the statement of Theorem V.3. Additionally, if these two equivalent statements are 
true, then V{9{t)) is non-increasing in A(7) for all 7 G [7min,7max], where 7min G [0,vr/2 — (fmaxl and 
7max £ ]^/2,vr] are given as unique solutions to inequality (30) with equality sign. Moreover, V{9{t)) is 
strictly decreasing in A (7) for all 7 G ]7min7 7max[- This concludes the proof of statement 1) and ensures 
that for every 0{O) G A(7max), there exists T > such that 0{t) G A{7r/2- ip^^^) for all t > T. Thus, the 
positive invariance assumption of Theorem V.l is satisfied, and statement 2) of Theorem V.3 follows. ■ 
In summary. Theorem V.3 presents sufficient conditions for the synchronization of the non-uniform 
Kuramoto model and is based on the bound (26). Condition (26) is a worst-case bound, both on the 
parameters and on the initial angles. In the remainder of this section, we aim at deriving a two-norm type 
bound and require only connectivity of the graph induced by P = and not necessarily completeness. 

^Besides tighter bounding of the right-hand side of (29), the proof can alternatively be continued by adding and subtracting the 
coupling with zero phase shifts in (29) or by noting that the right-hand side of (29) is a convex function of dk G [Oi, 6m] that 
achieves its maximum at the boundary 9^ £ {Oi, Om}- If the analysis is restricted to 7 G [0, 7r/2], the term bmk can be dropped. 
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We start our discussion with some preliminary notation and concepts. The following analysis is formally 
carried out for the complete graph, but, without loss of generality, we assume that some weights Pij = Pji 
can be zero and the non-zero weig hts P = P'^ induce a connected graph. Let H G M"("-i)/2xn jj^g 
incidence matrix of the complete graph with n nodes and recall that for a vector x G M" the vector of 
all difference variables is Hx = (x2 — xi, . . . ). The phase difference dynamics (24) (with the sinusoidal 
coupling expanded as in (18)) can be reformulated in a compact vector notation as 

j-^He = HD-^uj - HD-^H'^ diag(Pij cos((/Py)) sin{He) - HX , (31) 

where X G is the vector of lossy coupling with components Xi = Yl^=ii-^ij / ^i) cos{6i — 6j) 

and sin{H9) is the multivariable sine. The set of differential equations (31) is well defined on T": the 
left-hand side of (31) is the vector of frequency differences H6 = {02 — 9i, . . . ) taking values in the 
tangent space to T", and the right-hand side of (31) is a well-posed vector- valued function of 6* G T". 

With slight abuse of notation, we denote the two-norm of the vector of all geodesic distances by 
Il-f^^ll2 ~ (Si ~ ^iP)^^^> and aim at ultimately bounding the evolution of ||iJ0(t)||2. Following 

a classic Kuramoto analysis [33], [28], [27], [25], we note that the non-uniform Kuramoto model (8) 
with w = constitutes a Hamiltonian system with the Hamiltonian U{6)\uj=o defined in equation (6). 
An analysis of (31) by Hamiltonian arguments is possible, but results in very conservative conditions. 
In the recent Kuramoto literature [32], [33], a different Lyapunov function considered for the uniform 

1 1 2 

Kuramoto model (2) evolving on M" is simply H^f^Hg. Unfortunately, in the case of non-uniform rates 
Di this function's Lie derivative is sign-indefinite. However, it is possible to identify a similar Lyapunov 
function that has a Lie derivative with symmetric coupling. Consider the function W : T" — )• M defined by 

^(^) = \ EL E ■=! D^D, \e, - . (32) 

A Lyapunov analysis of system (31) via the Lyapunov function W leads to the following theorem. 

Theorem V.5 (Synchronization condition II) Consider the non-uniform Kuramoto model (8), where the 
graph induced by P = P^ is connected. Let H G ]K"("-i)/2xn incidence matrix of the complete 

graph and assume that the algebraic connectivity of the lossless coupling is larger than a critical value, i.e., 

\\HD-^uj\\^ + ^ YTj=i 5f sin(v?ij), . . . , Ej=i §: sin(v9„j) 



\2{L{Pij COs{ipij))) > \ 



critical 



cos((^ma)t) ('«/n)a/ maxj^j{A^j} 

(33) 



where k := X]fe=i <^nd a := ^minj^jlDjDj}/ maXj^jlDjDj}. 

Accordingly, define 7max G \'^/'^~ V^max,?!"] and 7min S [0,7r/2 — ipmax[ to be the unique solutions to 
the equations sinc(7max)/ sinc(7r/2 - c^max) = sin(7min)/ cos(99max) = Acriticai/A2(-^-(-Pji cos((^ij))). Then, 
1) phase cohesiveness: the set {9 G A(7r) : H^f^Hg < 7} is positively invariant for every 7 G 
[7min, 07max], and each trajectory starting in {6 G A(7r) : ||i70(O)||2 < a7max} reaches {9 G 
A(7r) : ||-ff^||2 < 7min}; and 
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2) frequency synchronization: /or every 9{0) G A(7r) with ||ii'6'(0)||2 < a7max the frequencies 6i{t) 
synchronize exponentially to some frequency O^o £ [^min(O), 0max(O)]. Moreover, if ifmax = 0, then 
6oo = ^ and the exponential synchronization rate is no worse than Afe as defined in equation (19). 

Remark V.6 (Physical interpretation of Theorem V.5:) In condition (33), ||[. . . , Yj^=i §^ sin(<^ij), • • -jjlg 
is the two-norm of the vector containing the lossy coupling, ||//Z)~^Lt;||2 = ||(w2/-D2 — '^iZ-Di) • • • )ll2 
corresponds to the non-uniformity in the natural frequencies, \2{L{Pij cos{ipij))) is the algebraic connec- 
tivity induced by the lossless coupling, cos((/?max) = sin(7r/2 — (pmax) reflects again the phase cohsiveness 
in A(7r/2 — v?max)> and {K/n)a/ maxj^jlDiDj} weights the non-uniformity in the time constants Di. 
The gap in condition (33) yields a again practical stability result determining the initial and ultimate 
phase cohesiveness. Condition (33) can be extended to non-reduced power network models [42]. □ 

Remark V.7 (Reduction of Theorem V.5 to classic Kuramoto oscillators:) For classic Kuramoto 
oscillators (2), condition (33) reduces to > identical '■— II ^'^112' which is a more conservative bound than 
K > identical = ^max " '^min presented in (27). It follows that the oscillators synchronize for ||ff0(O)||2 < 
7max and are ultimately phase cohesive in 11^^6*112 < 7min> where 7niax £ ]7r/2,7r] and 7min S [0,7r/2[ are 
the unique solutions to (vr/2) sinc(7max) = sin(7min) = -fir*riticai/-^- The Lyapunov function yV{6) reduces 
to the one used in [32], [33] and can also be used to prove [32, Theorem 4.2] and [33, Theorem 1]. □ 

Recall from Section I that angular differences are well defined for 6 G A(7r). Hence, for 9 G A(7r), 
the vector of phase differences is HO = {02 - 9i, ■■ . ) eM"^""^)/^, and the function W defined in (32) 
can be rewritten as the function H6 ^ W{H9) defined by 

^(^) = 7 . E" . ^^^i 1^^ - ^jl' = liHef <img{D,Dj){He) =: W{He) . (34) 
The derivative of W{H9) along trajectories of system (31) is then given by 

W{H9) = {H9f diag{DiDj)HD-^uj - {HOf diag{DiDj)HX 

- {H9)'^ diag{DiDj)HD-^H^ diag(Py cos(^y)) siii{H9) . (35) 
A component-wise analysis of the last term on the right-hand side of (35) yields a "diagonal" simplification. 

Lemma V.8 Let P = e W'", 9 e A(7r), and k := Y,k=i ^k- Then it holds that 
{H9f diag{DiDj)HD-^H^ diag(Pij cos((/7,j)) sin{HO) = K{H9f diag{Pij cos{ipij)) sin{H9). (36) 
Proof: The left-hand side of equation (36) reads component-wise as 
E, Efc(^^-^i)(^^'^' cos(^2fc)I?,) sin{9i-9k) + J] . J] J^2-%)(P,x. cos{ipjk)Di) sm{9k-9j) , 
where all indices satisfy /c G {1, . . . , n}. An manipulation of the indices in both sums yields 
E, Efe E ■i^^-^k){Pij cosiipij)Dk) sm{9i-9j) + E^ E • Y.S^'^-^^^^^'j cos(^y)Z)fc) sm{9i-9,) . 
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Finally, the two sums can be added and simplify to '^i'^f^'^jiPij cos{ipij)Dk){Oi — 6j)sm(9i — 6j), 
which equals the right-hand side of equation (36) written in components. ■ 
The following lemma will help us to upper-bound the derivative W{H6) by the algebraic connectivity. 

Lemma V.9 Consider a connected graph with n nodes induced by A = A^ ^W^^^ with incidence matrix 
B and Laplacian L{Aij). For any X £ M"", it holds that {Bx)"^ diag{Ai j){Bx) > {X2{L{Aij)) / n) ||-Bx||2 • 



Proof: Let H be the incidence matrix of the complete graph. The Laplacian of the complete graph 
with uniform weights is then given by (n • I„ — Inl^) = H^H, and the projection of x € M" on the 
subspace orthogonal to 1„ is x± = [In — = {l/n)H^ H x. Consider now the inequality 

(Bx)'^ diag{Aij)iBx) = x'^ B'^ diag{Aij)Bx = x^ L{Aij)x 

> A2(L(A.,)) llxJI^ = \\H^Hx\\l = ^^^^^^iHxfHH^iHx) . 

In order to continue, first note that HH^ and the complete graph's Laplacian H^H have the same 
eigenvalues, namely n and 0. Second, range(-ff) and ker(//^) are orthogonal complements. It follows 

rT-i / \ ii2 ii2 ii2 

that [Hx) HH (Hx) = n \\Hx\\2- Finally, note that ||-ffx||2 > ||^a:||2 ^rid the lemma follows. ■ 
Given Lemma V.8 and Lemma V.9 about the time derivative of W{H9), we are now in a position to 
prove Theorem V.5 via standard Lyapunov and ultimate boundedness arguments. 

Proof of Theorem V.5: Assume that 6^(0) G S{p) := {9 € A(7r) : \\H9\\2 < p} for some p £]0,tc[. In 
the following, we will show under which conditions and for which values of p the set S{p) is positively 
invariant. For 9 G S{p) and since < ||/?0||2, it follows that 9 G A(/j) and 1 > sinc(6'j — 9j) > 

sinc(p). Thus, for 9 G S{p), the inequality {9i — 9j) sm{9i — 9j) > {9i — 9j)'^smc{p) and Lemma V.8 
yield an upper bound on the right-hand side of (35): 

W{H9) < {H9fdiag{DiDj)HD-^u-{H9fdiag{DiDj)HX-Ksmc{p){H9fdiag{PijCos{ipij)){H9) . 



Note that ||//X||2 is lower bounded as 



--■.X. 



This lower bound X together with Lemma V.9 leads to the following upper bound on W{H9): 
W{H9) < \\H9\\2maiCi^j{DiDj}{\\HD-^u\\^ + X) - (k/u) smc{p)\2{L{Pij cos{ipij)) \\H . (37) 

Note that the right-hand side of (37) is strictly negative for 

_ maxi^,{D^Dj}i\\HD-^u;\\^ + X) 

In the following we apply standard Lyapunov and ISS arguments. Pick p £]0,p[. If 

^ max,^,{AI^,}(||ffI^-^a^||2 + ^) ^3^^ 
^ {K/n)smc{p)X2{L{Pij cos{(f ij))) 
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then for all ||-ff^||2 S [/^, /o], the right-hand side of (37) is upper-bounded by 

Wine) < -(1 - ific/f^)) ■ {K/n)smc{p)X2iLiPijCos{ipij))) \\He\\l . 

Note that W{H6) defined in (34) can easily be upper and lower bounded by constants multiplying ||i70||2: 

mmi^j{DiDj} \\H0\\l < 2 ■ W^HO) < maxi^j{D,Dj} \\He\\l . (39) 

To guarantee the ultimate boundedness of H9, two sublevel sets of W{H9) have to be fitted into {H9 : 
\\H9\\2 € where W{H9) is strictly negative. This is possible if [51, equation (4.41)] 



H < y miui^jlDiDj}/ maxi^j{DiDj} ■ p = ap . (40) 

Ultimate boundedness arguments [51, Theorem 4. 18] imply that, for every ||//^(0)||2 < q/3, there is T> 
such that \\H9{t)\\2 is strictly decreasing for t £ [0,T] and \\H9{t)\\2 < p/a for all t > T. If we choose 
p = with 7 g]0, 7r/2 — (/?max]> then equation (40) reduces to p > 7 and (38) reduces to the condition 

X2[L{Pij COs{ipij))) > . , . . l-r-^ = Acritical : ' (^1) 

a'y [K / n) smc[p) 7smc(p) 
where Acdticai is as defined in equation (33). Now, we perform a final analysis of the bound (41). The 
right-hand side of (41) is an increasing function of p and decreasing function of 7 that diverges to 00 
as p t vr or 7 I 0. Therefore, there exists some (p, 7) in the convex set A := {(p, 7) : p G]0, 7r[, 7 S 
]0, 7r/2— cpmax] , 7 < p} satisfying equation (41) if and only if equation (41) is true at p = 7 = 7r/2— (pmax, 
where the right-hand side of (41) achieves its infimum in A. The latter condition is equivalent to inequality 
(33). Additionally, if these two equivalent statements are true, then there exists an open set of points in A 
satisfying (41) , which is bounded by the unique curve that satisfies equation (41) with the equality sign, 
namely /(p,7) = 0, where / : A ^ M, /(p,7) := 7sinc(p)/ cos((pmax) - Acriticai/A2(-^^(^'ii cos((pij))). 
Consequently, for every (p, 7) G {{p,l) G A : f{p,"f) > 0}, it follows for 117/0(0)112 — '^P ^^^^ there 
is T > such that ||//6'(t)||2 < 7 for all t > T. The supremum value for p is obviously given by 
Pmax G ]7r/2 — ipmaxi 71"] solving the equation /(pmaxi tt/2 — (pmax) = and the corresponding infimum of 
7 by 7mm e [0, tt/2 - (pmax[ solving the equation /(7min, 7min) = 0. 

This proves statement 1) (where we replaced pmax by 7max) and shows that there is T > such that 
\\H9{t)\\^ < \\H9{t)\\2 < tt/2 - ip^^^ for all t > T. Statement 2) follows then from Theorem V.l. ■ 

C. Phase Synchronization 

For identical natural frequencies and zero phase shifts, the practical stability results in Theorem V.3 
and Theorem V.5 imply 7niin | 0, i.e., phase synchronization of the non-uniform Kuramoto oscillators (8). 

Theorem V.IO (Phase synchronization) Consider the non-uniform Kuramoto model (8), where the graph 
induced by P has a globally reachable node, (/Jmax = 0, and uJi/Di = u) for all i G {1, . . . , n}. 
Then for every 9{0) £ A(7) with 7 € [0, 7r[, 
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1) the phases 9i{t) synchronize exponentially to Ooo{t) S [0min(O), ^max(O)] + u)t, and 

2) if P = P^, then the phases 9i (t) synchronize exponentially to the weighted mean angle^ 9oo (t) = 

DiOiijS)/ Di + ujt at a rate no worse than 

Aps = -X2{L{Pj)) Sinc(7) COs(Z(Dl, l))VAnax . (42) 

The worst-case phase synchronization rate Aps can be interpreted similarly as the terms in (19), where 
sinc(7) corresponds to the initial phase cohesiveness in A(7). For classic Kuramoto oscillators (2) 
statements 1) and 2) can be reduced to the Kuramoto results found in [19] and Theorem 1 in [33]. 

Proof of Theorem V.IO: First we proof statement 1). Consider again the Lyapunov function V{9{t)) 
from the proof of Theorem V.3. The Dini derivative in the case (fmnK = and Wj/Dj = a) is simply 

D+V{9{t)) = - Y^l^^ sin(0^(t) - 9,{t)) + ^ sin{9,{t) - 9,{t))) . 

Both sinusoidal terms are positive for 9{t) G A(7), 7 G [0,7r[. Thus, V{9{t) is non-increasing, and A(7) 
is positively invariant. Therefore, the term aij{t) = {Pij / Di) smc{9i{t) — 9j{t)) is strictly positive for 
all t > 0, and after changing to a rotating frame (via the coordinate transformation 9 ^ 9 — ujt) the 
non-uniform Kuramoto model (8) can be written as the consensus time-varying consensus protocol 

^ii^) = - - ^j(O) , (43) 

Statement 1) follows directly along the lines of the proof of statement 1) in Theorem V.l. In the case of 
symmetric coupling P = P^ , the phase dynamics (43) can be reformulated as a symmetric time-varying 
consensus protocol with strictly positive weights Wij{t) = Pij sinc(6'j(t) —9j{t)) and multiple rates Di as 

^^D9 = -L{w,j{t))9, (44) 

Statement 2) now follows directly along the lines of the proof of statement 2) in Theorem V.l. ^ ■ 
The main result Theorem III.2 can be proved now as a corollary of Theorem V.3 and Theorem IV.2. 
Proof of Theorem 111.2: The assumptions of Theorem III.2 correspond exactly to the assumptions of 
Theorem V.3 and statements 1) and 2) of Theorem III.2 follow trivially from Theorem V.3. 

Since the non-uniform Kuramoto model synchronizes exponentially and achieves phase cohesiveness 
in A(7min) Si A(7r/2 — ipmax), it follows from Lemma IV.l that the grounded non-uniform Kuramoto 
dynamics (1 1) converge exponentially to a stable fixed point 5oo- Moreover, (5(0) = grnd(0(O)) is bounded 
and thus necessarily in a compact subset of the region of attraction of the fixed point 5oo. Thus, the 
assumptions of Theorem IV.2 are satisfied. Statements 3) and 4) of Theorem III.2 follow from Theorem 
IV.2, where we made the following changes: the approximation errors (14)-(15) are expressed as the 
approximation errors (10) in 6'-coordinates, we stated only the case e < e* and t > 4 > 0, we reformulated 
h(5{t)) = D^^Q{9{t)), and weakened the dependence of e on Q.^ to a dependence on 0(0). ■ 

^This weighted average of angles is geometrically well defined for 6^(0) G A(7r). 

''The proof of Theorem V.5 can be extended for HD~^lo — and X = to show statement 2) of Theorem V.IO with a 
slightly different worst-case synchronization frequency than (42). 
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VI. Simulation Results 

Figure 3 shows a simulation of the power network model (3) with n = 10 generators and the 
corresponding non-uniform Kuramoto model (8), where all initial angles ^(0) are clustered with exception 
of the first one (red curves) and the initial frequencies are chosen as ^(0) G uni(— 0.1, 0.1) rad/s, i.e., 
randomly from a uniform distribution over [—0.1, 0.1]. Additionally, at two-third of the simulation interval 
a transient high frequency disturbance is introduced at ujn~i (yellow curve). For illustration, relative 
angular coordinates are defined as 6i{t) = 6i{t) — 6n{t), i E {1, . . . , n— 1}. The system parameters satisfy 
e uni(0, 10), Pij G uni(0.7, 1.2), and tan((^ij) G uni(0, 0.25) matching data found in [11], [40], [41]. 

For the simulation in Figure 3(a), we chose Mj G uni(2, 12)s/(27r/o) and Di G uni(20, 30)s/(27r/o) 
resulting in the rather large perturbation parameter e = 0.58 s. The synchronization conditions of Theorem 
III.2 are satisfied, and the angles 5{t) of the non-uniform Kuramoto model synchronize very fast from the 
non-synchronized initial conditions (within 0.05 s), and the disturbance around t = 2s does not severely 
affect the synchronization dynamics. The same findings hold for the quasi-steady state h{5) depicting 
the frequencies of the non-uniform Kuramoto model, where the the disturbance at angle n — 1 (yellow 
curve) acts directly without being integrated. Since e is large the power network trajectories {5{t),6{t)) 
show the expected underdamped behavior and synchronize with second-order dynamics. As expected, the 
disturbance at t = 2 s does not affect the second order power network J-dynamics as much as the first- 
order non-uniform Kuramoto ^-dynamics. Nevertheless, after the initial and mid-simulation transients the 
singular perturbation errors 5{t) — 5{t) and 9{t) — h{5{t)) quickly become small and ultimately converge. 

Figure 3(b) shows the exact same simulation as in Figure 3(a), except that the simulation time is halved, 
the inertia are Mj G uni(2, 6)s/(27r/o), and the damping is chosen uniformly as = 30s/(27r/o), which 
gives the small perturbation parameter e = 0.18s. The resulting power network dynamics {5{t),6{t)) 
are strongly damped (note the different time scales), and the non-uniform Kuramoto dynamics 5{t) and 
the quasi-steady state h{5{t)) have slower time constants. As expected, the singular perturbation errors 
remain smaller during transients and converge faster than in the weakly damped case in Figure 3(a). 

VII. Conclusions 

This paper studied the synchronization and transient stability problem for a power network. A novel 
approach leads to purely algebraic conditions, under which a network-reduced power system model is 
transiently stable depending on network parameters and initial phase differences. Our technical approach 
is based on the assumption that each generator is highly overdamped due to local excitation control. The 
resulting singular perturbation analysis leads to the successful marriage of transient stability in power 
networks, Kuramoto oscillators, and consensus protocols. As a result, the transient stability analysis of a 
power network model reduces to the synchronization analysis of non-uniform Kuramoto oscillators. The 
study of generalized coupled oscillator models is an interesting mathematical problem in its own right 
and was tackled by combining and extending different techniques from all three mentioned areas. 
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(a) Weakly damped simulation with e = 0.58 s 



(b) Strongly damped simulation with e = 0.18 s 



Fig. 3. Simulation of the power network model (3) and the non-uniform Kuramoto model (8) 



The presented approach to synchronization in power networks offers easily checkable conditions and 
an entirely new perspective on the transient stability problem. The authors are aware that the derived 
conditions are not yet competitive with the sophisticated numerical algorithms developed by the power 
systems community. To render our results applicable to real power systems, tighter synchronization 
conditions have to be developed, the region of attraction has to be characterized more accurately, and 
more realistic power network models have to be considered. The authors' ongoing work addresses the 
last point and extends the presented analysis to structure-preserving power network models. 

Finally, the revealed relationship between power networks, Kuramoto oscillators, and consensus algo- 
rithms gives rise to various exciting research directions at the interface between these areas. 
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VIII. Appendix: Alternative Synchronization Conditions 

In this appendix we briefly comment on alternative bounding methods in the proof of Theorem V.3 and 
how they affect the triplet of the synchronization condition (9), the estimate for the region of attraction 
A(7max)> and the ultimate phase cohesive set A(7niin). We will state only the essential parts of the theorem 
statements and the corresponding proofs. 

A. Pairwise Bounding 

The proof of Theorem V.3 can be continued from equation (29) by bounding the right-hand side of 
equation (29) for each single pair {m, €\ rather than for all m, £ G {1, . . . , n}. Such a pairwise bounding 
results in n(n — l)/2 pairwise synchronization conditions and the worst multiplicative gap over all 
conditions determines the estimates for the region of attraction A(7niax) and the ultimate phase cohesive 
set A(7min). In short, tighter bounds are traded off for complexity. The resulting theorem statement is as 
follows. 

Theorem VIII.l (Synchronization condition I) Consider the non-uniform Kuramoto-model (8), where 
the graph induced by P = is complete. Assume that the minimal lossless coupling of any oscillator 
pair {m, €\ to the network is larger than a critical value, i.e., for every m,i £ {!,..., n}, m ^ I, 



^re{m,Q\{fe} [ Di J 

^critical 1 . [ _ ^ 



^ ("^ sin(v?mfc) + ^ sm{ipik)] j . (45) 
fc=i ^ ™ ^ ^ / 

Accordingly, define jmin £ [0,vr/2 — 93max[ and 7max G ]'^/2,7r] as unique solutions to the equations 

-icritk 



siii(7min) = sin(7max) = cos(v3max) maXm,e{T''^f''''^ /Tmi}. Then . . . 



Proof of Theorem VIII.l: . . . [see proof of Theorem V.3] ... 
In summary, D^V{9{t)) in (29) can be upper bounded by the simple expression 

D+V{9{t))< ^-^-Y: ^ min {a,,} sin(7) + V 5,, + V 6^, . 

It follows that V{6{t)) is non-increasing for all 6{t) £ A(7) and for all pairs {m,£} if 

Trrd sin(7) > cos((^^ax)r^f ^' , (46) 

where Tmi and T™"'^''' are defined in (45). The left-hand side of (46) is a strictly concave function 
of 7 G [0, vr]. Thus, there exists an open set of arc lengths 7 including 7* = it/2 — ifmax satisfying 
equation (46) if and only if equation (46) is true at 7* = 7r/2 — ipmax with the strict inequality sign, 
which corresponds to condition (45) in the statement of Theorem VIII.l. Additionally, if these two 
equivalent statements are true, then V{9{t)) is non-increasing in A (7) for all 7 G [7miin7max]> where 
7min S [0, tt/2 — (fmaxi and 7max S ]7r/2, vr] satisify inequaUty (46) for all pairs {m, £} . . . ■ 



June 28, 2011 



DRAFT 



33 



B. Pairwise Concavity-Based Bounding 

The bounding in the proof of Theorem V.3 can be further tightened by cocavity-based arguments. This 
bounding results in n{n — l)/2 pairwise synchronization conditions and n(n — 1) equations determining 
the estimates for the region of attraction A(7max) and the ultimate phase cohesive set A(7min)- Again, 
tighter bounds are traded off for increasing complexity. The resulting theorem statement is as follows. 

Theorem VIII.2 (Synchronization condition I) Consider the non-uniform Kuramoto-model (8), where 
the graph induced by P = is complete. Assume that the minimal lossless coupling of any oscillator 
pair {m, i} to the network is larger than a critical value, i.e., for every m,i £ {!,..., n}, m ^ I, 

• r r, 1 1 1^ cos(<^ifc + ^max) 1 > T^'f"' := ^ _ ^ + max J \ ■ (47) 

^ ie{r?i/}\{fc} [Vi J Dm P>£ ie{m,e} I ^ A I 



k: 



Accordingly, for every pair {m,l} define -y^f^ € [0,n/2 — ipmax[<3nd y^^^^ E]7r/2,7r] as unique solutions to 



E{ Pik ■ ( rat \ 1 ■ f ^ik ■ , mi , \\ -p. 

, iG{m,Q\{fc} J frt «e{m,£}\{fc} [ A J 



critical 

ml 1 



(48) 



iG{m,Q\{fc} \Di J t^l »e{m,n\{fc} Di 

and let y^a, := max„_<.{7^f } and 7max := minm^47™/^}. Then . . . 

Proof of Theorem VIII.2: . . . [see proof of Theorem V.3] ... 
Written out in components D^V{6{t)) (in the non-expanded form (8)) takes the form 

D+Vm) = ^ _ ^ _ V ^ smiOmit) - Okit) + 'fmk) + ^ sin(0,(t) - e,{t) - ■ (49) 

J-'m i^l J^m 

k=l 

In the following we abbreviate the summand on the right-hand side of (49) as fk{Qk{t)) '■= %^ sm{6m{t) — 

— ipik) and aim at a least conservative bounding of fk{0k{t)). There 
are different ways to continue from here - the following approach is based on concavity. 

Since fk{Ok) is the sum of two shifted concave sine functions of 9k G [9e, dm] (with same period and 
shifted by strictly less than it/2), fk{Ok) is again concave sine function of 9^ G [9e,9m] and necessarily 
achieves its minimum at the boundary 9^ G {9£,9m}- If argming^gjg fk{9k) = then 

fk{dk{t)) > /fc (7) := sin(7 + (pmk) - ^ sm{(pik) , 
Um Ui 

and otherwise for argming^,g[g^ g^j fk{9k) = Gm 

fkiOkit)) > fkil) ■= ^ sin(7 - m) + sin(v3„fc) , 

^^i ^m 

where fl and /| are functions from [0, vr] to M. In the following let /| : [0, vr] — )■ M be defined by 

fkil) ■= min I sin(7 + ifmk), ^ sin(7 - cpik) \ - ^ sm{(pik) . 

I Um P>i ) Ul 

Since fk{9{t)) > f^i'j) for all 9k{t) £ [9m{t), 9(,{t)], the derivative (49) is upper-bounded by 



Dm Dp 
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It follows that V{6{t)) is non-increasing for all 9{t) G A(7) and for all pairs {m,(.} if for all {m,* 

1 ^ l^m \ , I Pik ■ , N 

mm < sm(7 + sm(7 — 09,1.) > > + max < > sm 09,1.) ; 



(50) 

Note that the minimizing summand on the left-hand side of (50) is ^vi^i<^{ra,iY\{k}{Pik sin(7 — ^pik)/ Di} 
for 7 < 7r/2, minjg|„j £|'^|^|{Pjfc sin(7 + for 7 > 7r/2, and it achieves its maximum value 

'^™'i&{m,i}\{k}{Pik cos{ipik) / Di} for 7 = tt/2. In particular, there exists an open set of arc lengths 7 for 
including 7* = 7r/2 — ipmax for which V{6{t)) is non-increasing in A(7) if and only if inequality (50) is 
strictly satisfied for 7* = 7r/2 — (pmax- In this case, define j^^^ e [0, 7r/2 — ipmaxi and 7^^^ G ]7r/2, vr] as 
the two unique solutions to equation (50) with equality sign, which is equivalent to equation (48). Then 
V{9{t)) is non-increasing in A(7) for all 7 G hmiii^mlx]- Finally define 7min and 7max as the maximum 
and minimum values of 7,^^ and j^^^ over all pairs {m, £} . . . ■ 

C. Adding and Subtracting the Lossless Coupling 

The last possible bounding we explore is restricted to the set of initial conditions in A(7r/2 — 'Pmm) 
and results in a simple, scalar, and intuitive but very conservative synchronization condition. Instead of 
bounding the right-hand side of D~^V{6{t)) directly, we add and subtract the lossless coupling in the 
proof of Theorem V.3. The resulting theorem statement is as follows. 

Theorem VIII.3 (Synchronization condition I) Consider the non-uniform Kuramoto-model (8), where 
the graph induced by P = is complete. Assume that the minimal coupling is larger than a critical 
value, i.e., for every i,j G {1, . . . , n} 

-Pmin > -Pcriticai := — T (max ^ - ^ +maxV' ^sin((/7ij) ) . (51) 

Accordingly, define 7min = arcsin(cos(v3max)-Pcriticai/^'min) taking value in [0, 7r/2 - (/?max[ and 7max = 
7r/2 - (/9max- Then ... 

Proof of Theorem VIU.S: . . . see proof of Theorem V.3 . . . 
Written out in components (in the non-expanded form (8)) D^V{9{t)) takes the form 

D+V{e{t)) = ^ - ^ - V sinieUt) - Ok{t) + ^mk) - ^ sin(e,(t) - Ok{t) + m)l • 



(52) 

Adding and subtracting the coupling with zero phase shifts yields 

— sm(6'm(t) - 6'fc(t)) + — 

P-mk 



Dm Dp ^k=l V Dm Di 

sin(6'm(t) - dkit) + (frnk) - sm{9jn{t) - 9k{t)) 

+ yl .Tr( (t) - (t) ) + m ) + MOk it) - Oe (t) ) 



June 28, 2011 



DRAFT 



35 



Since both sinusoidal terms in the first sum are strictly positive, they can be lower-bounded as 

^ smiOmit) - Okit)) + ^ sin(0fc(t) - 0,(0) > ^ ( MOm{t) - 0fc(t)) + sin(0fc(t) - 0,(0)) • 

In the following we apply classic trigonometric arguments from the Kuramoto literature [32], [24], [39]. 
The identity sin(x) + sin(y) = 2 sin(^i^) cos(^^) leads to the further simplifications 

sin(0„(t) - e,{t)) + sin(0fc(t) - Oeit)) = 2 sin(^ ^"^'^ ~ ^^^'^ ) ^^^(^ O^it) + 9iit) _ ^^^^^ 

sin(0^(t) - 9k{t) + ^mk) + sin(0fc(t) - e^{t)) = 2 sin(^) cos(0„(t) - 0fc(t) + ^) , 

sin(0^(t) - Quit) + m) + sin(0fc(t) - Om = 2 sin(^) cos{e,{t) - Okit) + ^) . 

Note that the right-hand side of (52) is a convex function of 9^ G [Oe, Om] (see Proof of Theorem VIII.2) 
and accordingly achieves its maximum at the boundary for 6k G {Oi, 6m} Therefore, D~^V{9{t)) is upper 
bounded by 

\Di Dj i>max ^^==1 ^2/ V2/ 

E" -Pmfc„ . (^mk\ ( , V'mfcX , -Pft „ . ('^lk\ ( ^lk\ 
— — 2 sm cos 7 H + > ——2 sm cos . 

k=i Dm \ 2 J V ' 2 / ^k=i Di> \ 2 ) \ 2 ) 

Note that the second sum is strictly negative for 7 G [0,7r/2 — (/?max[ and can be neglected. Moreover, 

in the third sum the maximum over all nodes £ can be taken. Reversing the trigonometric identity from 

above as 2sin(x) cos(y) = sin(x — y) + sm{x + y) yields then the simple expression 



D+V{e{t)) < max ^ - ^ - Yll , ^111(7) + max^" ^ s\n{ipik) . 



It follows that the length of the arc formed by the angles is non-increasing in A (7) if 

Pmin Sin(7) > Pcritical COs{Lpm^y) , (53) 

where Pcriticai is as stated in equation (51). The left-hand side of (53) is a strictly increasing function of 7 G 
[0, 7r/2 — (/9max[- Therefore, there exists some 7* G [0, 7r/2 — i/9max[ satisfying equation (53) if and only if 
equation (53) at 7 = vr/2 — (/9max is true with the strict inequality sign, which corresponds to equation (51). 
Additionally, if these two equivalent statements are true, then there exists a unique j^m G [0, 7r/2 — (pmax[ 
that satisfies equation (53) with the equality sign, namely 7min = arcsin(cos((/?max)-Pcriticai/-fmin)- ■ ■ ■ 
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